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Abstract 

The problem of evaluating potential integrals on planar triangular el- 
ements has been addressed using a polar coordinate decomposition. The 
resulting formulae are general, exact, easily implemented, and have only 
one special case, that of a field point lying in the plane of the element. 
Results are presented for the evaluation of the potential and its gradients, 
where the integrals must be treated as principal values or finite parts, for 
elements with constant and linearly varying source terms. These results 
are tested by application to a single triangular element to the evaluation 
of the potential gradient outside the unit cube. In both cases, the method 
is shown to be accurate and convergent. 

1 INTRODUCTION 

A basic operation in any code for the Boundary Element Method (BEM) is 
the evaluation of potential integrals on elements, whether in the solution of the 
integral equation, or in the evaluation of field quantities. If we consider the 
Laplace problem, <j> the potential external to a surface S is given by the integral 
formulation: 



< / \ f \ <9G(x, Xl ) 



(1) 



where x indicates position, subscript 1 variables of integration on the surface S 
and n the outward pointing normal to the surface. The Green's function G is 
given by 

S(x;x l) = ^, (2) 

R = |x-xi|. 

Given the surface potential <j> and gradient d<j>/dn, the potential, and, after 
differentiation, the gradient(s), can be evaluated at any point in the field. Also, 
given a boundary condition for <j> and/or d(j>/dn on S, the integral equation can 
be solved for x € S. 

In any case, the method of solution remains the same: the surface S is 
divided into elements and suitable shape functions are used to interpolate the 



potential on these elements. The integral equation is transformed to a linear 
system in the element potentials, with the influence coefficients determined by 
the potential generated by each element at each node of the surface mesh. This 
leads to the requirement to evaluate integrals I and dl/dn where: 



where E is the surface of an element and (£, 77) is a coordinate system local to 



Numerous techniques have been proposed for the evaluation of / and its 
derivatives. Many of these have been numerical and have often focussed on the 
question of evaluating the integral when the evaluation point is close to, but 
not on, the element [1-3, for example]. For the Laplace equation, however, it is 
possible to exactly evaluate the potential integrals analytically and a number of 
techniques have been published using this approach [4-10, for example], which 
has the advantages of being exact and of being less prone to numerical errors 
introduced by singularities and near-singularities, and, often, more efficient than 
numerical quadrature. 

This paper introduces a method for the evaluation of potential integrals on 
flat triangular elements, in practice the most common problem encountered in 
BEM integration. The motivation for this work is to simplify the implementa- 
tion of the analytical formulae. The numerous results which have been published 
hereto are, obviously, algebraically equivalent in that they give, or should give, 
the same answer for any given combination of element and field point. In prac- 
tice, however, they are not numerically equivalent and, in particular, they have 
different special cases which must be handled. For example, the formulae of 
Newman [7] , based on the use of Green's theorem to reduce the surface integral 
to a sequence of line integrals on the element boundary, have a special case 
when the field point is collinear with an edge of the triangle. In Salvadori's 
work, there are a number of special cases which must be identified and handled 
separately depending on particular combinations of parameters [10, Figure 5]. 

In this paper, we present formulae based on a polar coordinate decomposition 
of the basic integrals which simplifies the integrands to the point where they 
can be easily evaluated using standard relations and tables and which require 
handling of only one special case, that when the field point lies in the plane 
of the element. The method can be generalized to source terms of any order, 
and to gradients of the potential, which give rise to strongly singular integrands 
which require treatment as hypcrsingular integrals. 

2 Potential integrals 

The approach adopted in developing formulae for the evaluation of potential 
integrals is similar to that used in previous work: the integral is evaluated 
on a reference triangle into which the real element can be transformed. In 
this case, the reference triangle lies in the plane z = 0, requiring a rotation 




(3) 
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Figure 1: Reference triangle for quadrature 



of the coordinate system, and has one vertex at the projection of the field 
point into that plane, requiring a decomposition of the original element into 
a set of subtriangles. We begin by developing formulae for the integrals over 
the reference triangle, before explaining how these formulae can be applied to 
general elements. 

Figure 3 shows the reference triangle and associated notation. The triangle 
lies in the plane z = 0, with one vertex at the origin. The two sides which meet 
at the origin have length n and r 2 and subtend an angle O. The third side is 
given by x = n + ay with: 



r 2 cos 9 — n 
r 2 sin 6 



(4) 



Integrations will be performed by transforming the double integrals in into inte- 
grals in the polar coordinate system (r, 9) centred on the origin. In the triangle 
coordinate system, the field point is placed at (0, 0, z). 

The integrals which are required arc of the following forms, with the factor 
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(5a) 
(5b) 
(5c) 
(5d) 
(5e) 
(5f) 



In these integrals, f{x\,y\) is a source term which varies over the element and 
subscript 1 denotes a variable of integration. The basic integral / is the integral 
over a general triangle in the plane z = of the Green's function weighted on 
the source term. Differentiation with respect to z, Equation 5b, corresponds to 
taking the normal derivative and derivatives with respect to (x, y, z) are used 
to find the gradient of potential, e.g. the velocity in potential flow problems. 
Further differentiation yields the hypersingular formulation used in certain ap- 
plications to avoid numerical difficulties introduced by spurious eigenvalues [11]. 

The integrals of Equations 5a can be evaluated on the reference triangle of 
Figure 3 as follows. Considering a general integral: 



= jj f{x 1 ,y l )R 1 dxidyi 



6 r r(0) 



Jo 



where 
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I can be evaluated analytically using tabulated relations. 

The first case considered is that of an clement with constant source distri- 
bution / (£1,2/1) = 1: 
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Under a change of variables: 
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(10) 
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which can be evaluated with a further change of variables from 6 to + <fi: 



7 + 2 
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where: 
A 2 



1 - a 2 sin 2 I 



fig 7 
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and R = /3 
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(12) 



rf + z 2 (l + a 2 
l + a 2 



Integrals of the form of Equation 12 are well understood and extensively tabu- 
lated [12, 13]. To simplify the statement of future results, we write the integral 
as: 



/ = 



1 



7 + 2 



P +2 (4^U^ e + <t>) - /f_ + 7 2 l 2 (a, 0)) - \zV +2 e] , (is) 



where 



I ( £ ) n (a,6)= sin m 6>'cos"0'A p d0'. 
Jo 



(14) 



This result is the general form for the integral of any integer power of R over 
the reference triangle. In particular, subject to triangle decomposition and 
rotation, it is the only integral required in a BEM which employs constant 
source elements, as it yields the integrals of G and dG/dn, and includes the 
integral of 1/R 3 considered by [9] . The only special case which need be handled 
is that for z = 0, when the field point lies in the element plane, which will be 
discussed in Section 2.1. 

Higher order elements can be handled in a similar fashion. In the case of 
elements with linear variation in f(xi,y\), the integrals to be evaluated are of 
the form: 
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In the case m = 0, n = 1, this yields: 
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upon integration by parts. Likewise, 
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d0, 
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(19) 

Finally higher order elements can be reduced to combinations of the integrals 
for constant and linear source triangles using the identity r 2 = R 2 — z 2 . For 
second order elements: 

r e+4> 
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which reduces to the standard integrals previously introduced: 
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Explicit expressions for the main results required for constant and linear ele- 

(v) 

ments are given in Table 1, with Im'n and J m , n given in the appendix, Tables 3 
and 4. 
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2.1 Field point in plane 

The formulae presented above are valid for all field points out of the triangle 
plane. The case of a point lying in the plane — the only special case which 
arises — must be handled separately for two reasons. The first is that the for- 
mulae presented so far break down numerically unless certain limits are taken 
explicitly; the second is that a number of the integrals in question are truly 
singular and must be interpreted as principal value or hypersingular integrals. 
It is more convenient to handle all of the in-plane cases together using the ap- 
proach of Brandao [14] who gives a convenient analysis for singular integrals 
of general form. For example, taking the integral of 1/-R 3 , and noting that for 
z = 0, r = R: 
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R 3 ' ~./n Vn r 



i r r i 

rdrd9= I I ^drd9. (23) 
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Interpreting the integral in r as a finite-part integral, using Brandao's method, 
it can be written: 
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The result of Brandao's which is applicable is: 
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so that, upon integration by parts: 
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Figure 2: Integration over a general triangle (left) by subdivision into three 
triangles centred at the origin (right). The triangle shown dashed on the right 
has negative orientation and its contribution is subtracted from that of the other 
two. 



Finally, for integrals containing terms R , the corresponding results are: 
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2.2 Integration over general triangles 

The preceding sections give exact results for integration over the reference tri- 
angle of Figure 3. Any triangular element can be transformed to a combination 
of triangles of reference shape by the following procedure. It is assumed that 
the coordinate system has already been transformed so that the triangle lies in 
the plane z = 0, a routine operation in evaluation of potential integrals. 

Figure 2 shows a general triangle 123 and the projection of a field point 
into the plane of the triangle at point 0. The system is decomposed into three 
subtriangles 012, 023, and 031 which are all of reference form. The integral 
over the triangle is then the sum of the integrals over the subtriangles. In 
summing the subtriangle integrals, account must be taken of the orientation of 
the triangles. For example, in Figure 2, the contributions of triangles 012 and 
023 are added, while that of 031 is subtracted. The simplest way to incorporate 
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Table 1: Summary of main results for constant and linear triangular elements 
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the orientations is to apply it to the calculation of the angle in Figure 3: 

-i (x 2 - x ).(xi - x ) 



= ± cos 



|x 2 - xo xi - x 



(33) 



where the sign is chosen to agree with the orientation (signed area) of the 
triangle. Orientations are computed using Shewchuk's robust adaptive predi- 
cates [15]. 

For constant source elements, the subtriangle integrals can be summed di- 
rectly. When there is a variation in the source term, two rotations must be 
applied, one by angle fa Equation 8, and one by angle tp: 



tp = tan 



-i yi - yo 



(34) 



x\ - x Q 

where subscript refers to the field point and subscript 1 to the first vertex of 
the triangle. Assembling these rotations gives, for the linear source: 
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where subscript i = 1, 2,3 refers to the subtriangles into which the element is 
decomposed. 



2.3 Summary of algorithm 

In summary, given a triangle lying in the plane z — and a field point at some 
general position, the potential integrals can be computed as follows: 

1. check if the field point lies in the plane z = 0; 

2. decompose the triangle into subtriangles as in Figure 2; 

3. for each triangle i — 1,2,3: 

(a) compute the orientation (signed area) of the triangle; 

(b) if the orientation is not zero, calculate a, ri, r 2 , 0i, and fa using 
Equations 4 and 8, taking account of the orientation; 

(c) evaluate the required integrals using the appropriate formulae of Ta- 
ble 1, depending on whether or not z = 0; 

(d) if required, apply the rotations of Equation 35. 

4. sum the results from the three subtriangles. 

If necessary, e.g. for the computation of gradients, a further rotation can then 
be applied to return to the global coordinate system. 
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(0,0.4) 



1 




(-0.5,-0.1) 



4 



(0.5,-0.1) 



Figure 3: Reference triangle for test of integration procedures. Points 1-5 are the 
projections onto the triangle plane of the field points. Points are 1: (—0.5, —0.1); 
2: (0.1,0.1); 3: (0.1, -0.101); 4: (-0.1, -0.099); 5: (0.3,0.2). 

3 Numerical tests 

Two main sets of tests have been performed to check the formulae presented. 
The first consists of the evaluation of sample integrals on a single triangle, with 
field points chosen to be representative of the cases likely to cause numerical 
difficulties. The second is a test of the estimation of the gradient of a known 
potential field outside the unit cube. This tests the method in a BEM code 
and assesses its ability to accurately integrate the strongly singular integrands 
which arise in computing gradients. 

3.1 Integrals on a sample triangle 

The triangular element for the first test is shown in Figure 3. The test integrals 
were evaluated at a range of values of z, over five different positions in the 
triangle plane, shown in Figure 3. For I\ and I2, < z < 1, while for I 3 , 
1/8 < z < 1, since the numerical integration routine will not give correct results 
at z = 0. These points were chosen to test for the cases most likely to be 
encountered in practice. Point 1 is directly over a triangle vertex; point 2 lies 
inside the triangle boundary; points 3 and 4 lie near the triangle boundary but 
just inside and outside it respectively; point 5 is well separated from the element. 
Three basic integrals were evaluated: 



where L^, 77), i = 1,2,3 are the linear shape functions for the triangle: 




(36b) 



(36a) 



(36c) 
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Table 2: Errors in evaluation of integrals for triangle of Figure 3 
Point £i £2 £3 



1 


1.5 x 10- 


-14 


1.4 X 10" 


-12 


4.4 x 10- 


-10 


2 


2.5 x 10" 


-li 


1.4 x 10" 


-09 


4.3 x 10" 


-07 


3 


4.9 x 10" 


-04 


6.9 x 10 


-09 


2.4 x 10" 


06 


4 


1.7 x 10" 


-05 


6.8 x 10 


-09 


2.4 x 10" 


06 


5 


4.6 x 10" 


-14 


2.9 x 10 


-12 


8.9 x 10" 


-10 




Figure 4: Unit cube and gradient test surface 



For comparison with a numerical method, the integrals were also evaluated 
using the approach of Hayami [3]. The error e,-, j — 1, 2, 3, is calculated as the 
maximum of the differences between the numerical and analytical values for Ij , 
for each value of i, and is shown in Table 2. 

From Table 2, the accuracy of the method is clearly displayed: the errors 
are comparable to machine precision, with the exception of the hypersingular 
case for points 2-4, where the numerical integration routine would be expected 
to break down, and for points 3 and 4 for I\, where the numerical method has 
difficulty in dealing with the very slender triangle between the evaluation point 
and the lower vertices of the triangle. 

3.2 Potential gradient calculation 

The second test was an evaluation of the gradient of a potential field outside 
the unit cube. The evaluation method was implemented in the free BEM code 
BEM3D [16], based on the GTS triangulated surface library [17], and used to 
find the gradient of the potential outside a unit cube centred on the origin, 
Figure 4. The boundary condition for potential and potential gradient was 
imposed using a point source placed inside the cube at (1/4, 1/4, 1/4) and the 
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-6 -5 -4 -3 -2 
\og 2 h 



Figure 5: Error versus discretization length for gradient of Laplace potential 
outside the unit cube. Symbols: computed error; solid line: e = 2.72ft 212 . 

field was evaluated as the gradient of Equation 1. The measure of error was the 
maximum difference in any of the three components of the gradient on a 32 x 32 
grid -2 < x < 2, -2 < y < 2, z = 0: 

/max|V0|. (38) 

The cube mesh was generated using GMSH [18], with the discretization length h 
varied as a parameter. 

The resulting error is shown in Figure 5 and demonstrates second order con- 
vergence with h, as might be expected for linear elements. Again, the accuracy 
of the method is confirmed. 

4 Conclusions 

A method and explicit formulae for the exact integration of potential integrals 
on planar triangular elements in the boundary element method has been pre- 
sented. The formulae are easily implemented using standard operations and 
testing against numerical quadrature has shown them to be accurate and, when 
implemented in a BEM code, convergent. 

A Basic integrals 

Expressions for Im]n and J m ,n can be found in standard tables [12, 13] and are 

( — 3) 

listed in Tables 3 and 4. Also, where Im,n is not given explicitly, it can be 
expressed in terms of lower order integrals [12, 2.581.2]: 

T<-V(n fi\ - sm^flcos"-^ m-1 n - 1 ( _i) 

(39) 



e = max 
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exact 
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Table 3: I&> 


(ft), constant of integration omitted 
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Table 4: J mn (0), constant of integration omitted 

Til Ti J mn 
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